import pandas as pd
import numpy as np
import geopandas as gpd # Vector geospatial operations
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
import datetime
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow. warnings.warn(
['lds-nz-road-centrelines-topo-150k-FGDB.zip', 'subnational-population-projections-2018base-2048.xlsx', 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip', 'statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip', 'statsnzregional-council-2021-clipped-generalised-FGDB.zip', 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']
First, read regional council bounds. This geometry will be used to clip NZ-wide datasets to just the region of interest, Auckland
%%time
REGC = gpd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)
epsg:2193 CPU times: user 1.59 s, sys: 92.1 ms, total: 1.68 s Wall time: 20.2 s
Load the LRIS Land Cover Database (downloaded in GDB format from https://lris.scinfo.org.nz/layer/104400-lcdb-v50-land-cover-database-version-50-mainland-new-zealand/)
%%time
df = gpd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")
CPU times: user 1min 28s, sys: 2.5 s, total: 1min 30s Wall time: 1min 30s
print(df.columns)
print(df.crs)
display(df.sample(5))
Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
dtype='object')
epsg:2193
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 231761 | Tall Tussock Grassland | Tall Tussock Grassland | Tall Tussock Grassland | Tall Tussock Grassland | Tall Tussock Grassland | 43 | 43 | 43 | 43 | 43 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30T00:00:00 | lcdb1000062986 | MULTIPOLYGON (((1827978.878 5664572.742, 18279... |
| 286399 | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | 54 | 54 | 54 | 54 | 54 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2015-06-30T00:00:00 | lcdb1000453515 | MULTIPOLYGON (((1781263.403 5625969.086, 17812... |
| 219917 | Sub Alpine Shrubland | Sub Alpine Shrubland | Sub Alpine Shrubland | Sub Alpine Shrubland | Sub Alpine Shrubland | 55 | 55 | 55 | 55 | 55 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb2000157563 | MULTIPOLYGON (((1611649.433 5390608.016, 16116... |
| 271790 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000205695 | MULTIPOLYGON (((1937691.648 5662380.164, 19377... |
| 333692 | Low Producing Grassland | Low Producing Grassland | Low Producing Grassland | Low Producing Grassland | Low Producing Grassland | 41 | 41 | 41 | 41 | 41 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30T00:00:00 | lcdb2000068052 | MULTIPOLYGON (((1236430.843 5030840.782, 12364... |
5 rows × 24 columns
%%time
df = gpd.clip(df, AKL)
CPU times: user 50.5 s, sys: 8.88 ms, total: 50.5 s Wall time: 50.5 s
df.sample(5)
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 345072 | Gorse and/or Broom | Gorse and/or Broom | Gorse and/or Broom | High Producing Exotic Grassland | High Producing Exotic Grassland | 51 | 51 | 51 | 40 | 40 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000412356 | POLYGON ((1779084.476 5881655.907, 1779073.885... |
| 1322 | Mangrove | Mangrove | Mangrove | Mangrove | Mangrove | 70 | 70 | 70 | 70 | 70 | ... | yes | no | no | no | no | no | Landcare Research | 2004-06-30T00:00:00 | lcdb1000182100 | POLYGON ((1754147.706 5940019.597, 1754156.167... |
| 478659 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000420376 | POLYGON ((1736616.736 5920808.002, 1736618.854... |
| 501381 | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | 56 | 56 | 56 | 56 | 56 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000129574 | POLYGON ((1732504.123 5936906.171, 1732487.177... |
| 393102 | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | 52 | 52 | 52 | 52 | 52 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000095134 | POLYGON ((1727359.405 5960844.635, 1727352.011... |
5 rows × 24 columns
df.Name_2018.value_counts()
Exotic Forest 3981 Indigenous Forest 3673 Manuka and/or Kanuka 2282 Broadleaved Indigenous Hardwoods 1788 Built-up Area (settlement) 1350 High Producing Exotic Grassland 1326 Mangrove 1151 Urban Parkland/Open Space 1099 Estuarine Open Water 441 Orchard, Vineyard or Other Perennial Crop 436 Short-rotation Cropland 362 Lake or Pond 326 Herbaceous Saline Vegetation 303 Low Producing Grassland 291 Gorse and/or Broom 287 Forest - Harvested 266 Sand or Gravel 252 Deciduous Hardwoods 201 Surface Mine or Dump 132 Mixed Exotic Shrubland 120 Herbaceous Freshwater Vegetation 118 Transport Infrastructure 107 River 15 Gravel or Rock 9 Flaxland 9 Fernland 2 Matagouri or Grey Scrub 1 Name: Name_2018, dtype: int64
These classes are far too detailed - simplify to just Urban, Vegetation, Water, Other
urban = [1,2,5]
vegetation = [68,69,71]
water = [0,20,21,22,45,46]
def simplify_classes(code):
if code in urban:
return 1, "Urban"
elif code in vegetation:
return 2, "Vegetation"
elif code in water:
return 3, "Water"
else:
return 4, "Other"
summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
print(year)
name_year = f"Name_{year}"
class_year = f"Class_{year}"
df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
display(df[[name_year, class_year]][df[class_year].isin(urban)].value_counts())
display(df[[name_year, class_year]][df[class_year].isin(vegetation)].value_counts())
display(df[[name_year, class_year]][df[class_year].isin(water)].value_counts())
summary.append(df[class_year + "_simplified_name"].value_counts())
1996
Name_1996 Class_1996 Urban Parkland/Open Space 2 1248 Built-up Area (settlement) 1 489 Transport Infrastructure 5 67 dtype: int64
Name_1996 Class_1996 Exotic Forest 71 3955 Indigenous Forest 69 3737 Deciduous Hardwoods 68 202 dtype: int64
Name_1996 Class_1996 Estuarine Open Water 22 448 Herbaceous Saline Vegetation 46 313 Lake or Pond 20 285 Herbaceous Freshwater Vegetation 45 117 River 21 15 Not land 0 3 dtype: int64
2001
Name_2001 Class_2001 Urban Parkland/Open Space 2 1167 Built-up Area (settlement) 1 748 Transport Infrastructure 5 74 dtype: int64
Name_2001 Class_2001 Exotic Forest 71 4495 Indigenous Forest 69 3717 Deciduous Hardwoods 68 206 dtype: int64
Name_2001 Class_2001 Estuarine Open Water 22 446 Herbaceous Saline Vegetation 46 312 Lake or Pond 20 305 Herbaceous Freshwater Vegetation 45 116 River 21 15 Not land 0 3 dtype: int64
2008
Name_2008 Class_2008 Built-up Area (settlement) 1 1082 Urban Parkland/Open Space 2 1080 Transport Infrastructure 5 108 dtype: int64
Name_2008 Class_2008 Exotic Forest 71 4562 Indigenous Forest 69 3695 Deciduous Hardwoods 68 202 dtype: int64
Name_2008 Class_2008 Estuarine Open Water 22 441 Lake or Pond 20 311 Herbaceous Saline Vegetation 46 310 Herbaceous Freshwater Vegetation 45 118 River 21 15 Not land 0 3 dtype: int64
2012
Name_2012 Class_2012 Built-up Area (settlement) 1 1136 Urban Parkland/Open Space 2 1091 Transport Infrastructure 5 105 dtype: int64
Name_2012 Class_2012 Exotic Forest 71 4345 Indigenous Forest 69 3690 Deciduous Hardwoods 68 201 dtype: int64
Name_2012 Class_2012 Estuarine Open Water 22 441 Lake or Pond 20 321 Herbaceous Saline Vegetation 46 306 Herbaceous Freshwater Vegetation 45 119 River 21 15 Not land 0 3 dtype: int64
2018
Name_2018 Class_2018 Built-up Area (settlement) 1 1350 Urban Parkland/Open Space 2 1099 Transport Infrastructure 5 107 dtype: int64
Name_2018 Class_2018 Exotic Forest 71 3981 Indigenous Forest 69 3673 Deciduous Hardwoods 68 201 dtype: int64
Name_2018 Class_2018 Estuarine Open Water 22 441 Lake or Pond 20 326 Herbaceous Saline Vegetation 46 303 Herbaceous Freshwater Vegetation 45 118 River 21 15 dtype: int64
pd.DataFrame(summary).plot.area()
<AxesSubplot:>
%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
ims = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
ax.set_title(year)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("land_use.gif", ims, fps=1)
with open('land_use.gif','rb') as file:
display(Image(file.read()))
cols = [f"Class_{year}_simplified_code" for year in years]
cols
['Class_1996_simplified_code', 'Class_2001_simplified_code', 'Class_2008_simplified_code', 'Class_2012_simplified_code', 'Class_2018_simplified_code']
%%time
geocube = make_geocube(
vector_data=df,
output_crs="epsg:2193",
measurements=cols,
resolution=(-100, 100),
fill=0,
)
geocube
CPU times: user 15.8 s, sys: 19.5 ms, total: 15.9 s Wall time: 15.9 s
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
spatial_ref int64 0
Data variables:
Class_1996_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2001_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2008_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2012_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2018_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])# Replace offshore areas (0) with 3 (Water)
geocube = geocube.where(geocube != 0, 3)
geocube.Class_2018_simplified_code.plot()
<matplotlib.collections.QuadMesh at 0x7fb9bc81a940>
for year in years:
print(year)
outfile = f"output/land_use_{year}.tif"
if not os.path.isfile(outfile):
geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize
1996 2001 2008 2012 2018
%%time
TALB = gpd.read_file("input/statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip!territorial-authority-local-board-2021-clipped-generalised.gdb")
TALB = gpd.clip(TALB, AKL)
TALB.head()
CPU times: user 1.37 s, sys: 9.55 ms, total: 1.38 s Wall time: 1.42 s
| TALB2021_V1_00 | TALB2021_V1_00_NAME | TALB2021_V1_00_NAME_ASCII | LAND_AREA_SQ_KM | AREA_SQ_KM | Shape_Length | geometry | |
|---|---|---|---|---|---|---|---|
| 2 | 00300 | Kaipara District | Kaipara District | 3108.960347 | 3108.960347 | 934117.838186 | MULTIPOLYGON (((1718576.593 5980699.195, 17185... |
| 4 | 01200 | Hauraki District | Hauraki District | 1270.133527 | 1270.133527 | 311282.487609 | POLYGON ((1801726.041 5899678.947, 1803788.474... |
| 5 | 01300 | Waikato District | Waikato District | 4403.797301 | 4451.122298 | 680627.805117 | POLYGON ((1747419.052 5870470.171, 1747413.783... |
| 44 | 07601 | Rodney Local Board Area | Rodney Local Board Area | 2275.114945 | 2275.114945 | 920471.582397 | MULTIPOLYGON (((1717309.547 5972989.294, 17173... |
| 45 | 07602 | Hibiscus and Bays Local Board Area | Hibiscus and Bays Local Board Area | 110.097301 | 110.097301 | 162589.255306 | MULTIPOLYGON (((1751233.745 5948474.046, 17512... |
# From https://www.stats.govt.nz/information-releases/subnational-population-projections-2018base2048#map
pop = pd.concat(pd.read_excel(
"input/subnational-population-projections-2018base-2048.xlsx",
sheet_name=["Table 5", "Table 6"],
skiprows=5,
usecols="A,B,G",
names=["area", "year", "population"],
nrows=231
)).fillna(method='ffill').reset_index()
pop.replace("Maungakiekie-Tamaki local board area", "Maungakiekie-Tāmaki local board area", inplace=True)
pop.year = pop.year.astype(str)
pop
| level_0 | level_1 | area | year | population | |
|---|---|---|---|---|---|
| 0 | Table 5 | 0 | Far North district | 1996 | 54500 |
| 1 | Table 5 | 1 | Far North district | 2001 | 56400 |
| 2 | Table 5 | 2 | Far North district | 2006 | 57500 |
| 3 | Table 5 | 3 | Far North district | 2013 | 60600 |
| 4 | Table 5 | 4 | Far North district | 2018 | 67900 |
| ... | ... | ... | ... | ... | ... |
| 457 | Table 6 | 226 | Franklin local board area | 2028 | 96900 |
| 458 | Table 6 | 227 | Franklin local board area | 2033 | 109400 |
| 459 | Table 6 | 228 | Franklin local board area | 2038 | 122000 |
| 460 | Table 6 | 229 | Franklin local board area | 2043 | 134500 |
| 461 | Table 6 | 230 | Franklin local board area | 2048 | 146900 |
462 rows × 5 columns
pop[pop.area.str.contains("local")].pivot(index='year', columns='area', values='population').plot()
<AxesSubplot:xlabel='year'>
pop = pop.pivot(index='area', columns='year', values='population')
pop
| year | 1996 | 2001 | 2006 | 2013 | 2018 | 2023 | 2028 | 2033 | 2038 | 2043 | 2048 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| area | |||||||||||
| Albert-Eden local board area | 86100 | 90000 | 96200 | 100000 | 103700 | 106900 | 111200 | 116500 | 121600 | 126300 | 130600 |
| Auckland | 1115800 | 1218300 | 1373000 | 1493200 | 1654800 | 1778700 | 1891800 | 2001800 | 2107000 | 2207800 | 2302900 |
| Devonport-Takapuna local board area | 51600 | 52300 | 55700 | 58500 | 60500 | 62700 | 65700 | 68800 | 71800 | 74800 | 77600 |
| Far North district | 54500 | 56400 | 57500 | 60600 | 67900 | 72400 | 75000 | 77300 | 79100 | 80600 | 81700 |
| Franklin local board area | 48400 | 52900 | 60600 | 68300 | 77700 | 85700 | 96900 | 109400 | 122000 | 134500 | 146900 |
| Gisborne district | 47200 | 45500 | 45900 | 47000 | 49500 | 51900 | 53100 | 54000 | 54600 | 55000 | 55200 |
| Great Barrier local board area | 1230 | 1100 | 940 | 950 | 960 | 1030 | 1080 | 1120 | 1150 | 1180 | 1200 |
| Hamilton city | 113500 | 121200 | 134800 | 150200 | 168600 | 183000 | 194400 | 205400 | 216000 | 226500 | 236600 |
| Hauraki district | 18550 | 18000 | 18300 | 18600 | 20700 | 21800 | 22200 | 22400 | 22300 | 22100 | 21800 |
| Henderson-Massey local board area | 82200 | 90900 | 103700 | 113500 | 124600 | 134100 | 142500 | 150500 | 158200 | 165700 | 172900 |
| Hibiscus and Bays local board area | 68000 | 75100 | 85200 | 94000 | 108500 | 117700 | 123500 | 126500 | 128900 | 130900 | 132400 |
| Howick local board area | 80900 | 97300 | 119100 | 135000 | 149400 | 160500 | 167800 | 174400 | 180400 | 185900 | 190900 |
| Kaipara district | 17800 | 17950 | 18550 | 20500 | 23700 | 25900 | 27200 | 28300 | 29100 | 29800 | 30300 |
| Kaipātiki local board area\n | 74700 | 78200 | 83600 | 87000 | 92900 | 97300 | 99400 | 101100 | 102500 | 103600 | 104700 |
| Kawerau district | 8120 | 7290 | 7150 | 6650 | 7460 | 7910 | 8000 | 8020 | 7970 | 7860 | 7720 |
| Manurewa local board area | 60100 | 69100 | 81600 | 87000 | 100900 | 109600 | 114000 | 118000 | 121700 | 124900 | 127600 |
| Matamata-Piako district | 30300 | 30300 | 31200 | 32900 | 35300 | 37000 | 37900 | 38700 | 39200 | 39500 | 39600 |
| Maungakiekie-Tāmaki local board area | 61700 | 66400 | 70400 | 73700 | 80500 | 86600 | 93800 | 101300 | 108800 | 116100 | 123000 |
| Māngere-Ōtāhuhu local board area\n | 61400 | 64600 | 72500 | 75300 | 82700 | 88700 | 92500 | 96300 | 99800 | 103100 | 105900 |
| Papakura local board area | 38000 | 38900 | 43100 | 48200 | 61100 | 70000 | 77000 | 81700 | 86100 | 90500 | 95000 |
| Puketāpapa local board area\n | 46100 | 49800 | 53900 | 56300 | 60900 | 64600 | 68800 | 73000 | 77100 | 81000 | 84600 |
| Rodney local board area | 39000 | 44100 | 51000 | 57300 | 69100 | 77600 | 87600 | 99600 | 111700 | 123800 | 135800 |
| Rotorua district | 66600 | 66900 | 68100 | 68400 | 74800 | 78900 | 80700 | 82200 | 83400 | 84200 | 84800 |
| South Waikato district | 25800 | 24200 | 23200 | 23200 | 24800 | 25800 | 26400 | 26800 | 27000 | 27100 | 27100 |
| Taupō district\n | 31600 | 32500 | 33400 | 34800 | 38600 | 40900 | 42000 | 42800 | 43300 | 43600 | 43800 |
| Tauranga city | 79900 | 93600 | 107000 | 119900 | 142100 | 156900 | 166300 | 175000 | 183300 | 191400 | 199100 |
| Thames-Coromandel district | 25400 | 25800 | 26700 | 27300 | 30700 | 32400 | 33100 | 33500 | 33500 | 33300 | 32800 |
| Upper Harbour local board area | 24300 | 33700 | 45100 | 56800 | 66800 | 75900 | 86000 | 95900 | 105800 | 115400 | 124900 |
| Waiheke local board area | 6680 | 7560 | 8150 | 8630 | 9360 | 9830 | 10250 | 10650 | 10900 | 11100 | 11200 |
| Waikato district | 52000 | 53700 | 59500 | 66500 | 78200 | 86100 | 92800 | 99300 | 105700 | 111900 | 117700 |
| Waipa district | 38400 | 40000 | 43700 | 48700 | 55000 | 59300 | 62100 | 64600 | 66900 | 68900 | 70700 |
| Waitematā local board area\n | 44900 | 52500 | 66600 | 81300 | 88500 | 94700 | 102800 | 111900 | 120700 | 129200 | 137400 |
| Waitomo district | 10000 | 9780 | 9680 | 9340 | 9630 | 9740 | 9730 | 9670 | 9540 | 9330 | 9070 |
| Waitākere Ranges local board area\n | 42100 | 44200 | 47500 | 50700 | 54200 | 56900 | 58300 | 59500 | 60300 | 60900 | 61100 |
| Western Bay of Plenty district | 35600 | 38900 | 42900 | 45400 | 53300 | 58100 | 60900 | 63300 | 65200 | 66700 | 68000 |
| Whakatāne district\n | 34200 | 34100 | 34500 | 34200 | 37100 | 38800 | 39300 | 39500 | 39500 | 39300 | 38900 |
| Whangārei district\n | 68400 | 70000 | 76500 | 83700 | 94100 | 101000 | 105500 | 109500 | 113100 | 116400 | 119300 |
| Whau local board area | 59700 | 66200 | 73100 | 76700 | 84100 | 89500 | 95000 | 100500 | 105900 | 111000 | 115900 |
| Ōpōtiki district\n | 9620 | 9490 | 9200 | 8780 | 9670 | 10250 | 10350 | 10400 | 10300 | 10150 | 9910 |
| Ōrākei local board area\n | 70900 | 73000 | 78400 | 83700 | 87700 | 91300 | 96200 | 101100 | 105800 | 110300 | 114500 |
| Ōtara-Papatoetoe local board area\n | 67800 | 70400 | 76600 | 80300 | 90500 | 97700 | 101500 | 104000 | 106000 | 107500 | 108600 |
| Ōtorohanga district\n | 9960 | 9590 | 9310 | 9590 | 10500 | 11050 | 11300 | 11550 | 11750 | 11900 | 12000 |
cols = pop.columns.tolist()
cols
['1996', '2001', '2006', '2013', '2018', '2023', '2028', '2033', '2038', '2043', '2048']
pop = pd.merge(TALB, pop, left_on=TALB.TALB2021_V1_00_NAME.str.lower().str.strip(), right_on=pop.index.str.lower().str.strip(), how="left")
pop[["TALB2021_V1_00_NAME"] + cols]
| TALB2021_V1_00_NAME | 1996 | 2001 | 2006 | 2013 | 2018 | 2023 | 2028 | 2033 | 2038 | 2043 | 2048 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Kaipara District | 17800 | 17950 | 18550 | 20500 | 23700 | 25900 | 27200 | 28300 | 29100 | 29800 | 30300 |
| 1 | Hauraki District | 18550 | 18000 | 18300 | 18600 | 20700 | 21800 | 22200 | 22400 | 22300 | 22100 | 21800 |
| 2 | Waikato District | 52000 | 53700 | 59500 | 66500 | 78200 | 86100 | 92800 | 99300 | 105700 | 111900 | 117700 |
| 3 | Rodney Local Board Area | 39000 | 44100 | 51000 | 57300 | 69100 | 77600 | 87600 | 99600 | 111700 | 123800 | 135800 |
| 4 | Hibiscus and Bays Local Board Area | 68000 | 75100 | 85200 | 94000 | 108500 | 117700 | 123500 | 126500 | 128900 | 130900 | 132400 |
| 5 | Upper Harbour Local Board Area | 24300 | 33700 | 45100 | 56800 | 66800 | 75900 | 86000 | 95900 | 105800 | 115400 | 124900 |
| 6 | Kaipātiki Local Board Area | 74700 | 78200 | 83600 | 87000 | 92900 | 97300 | 99400 | 101100 | 102500 | 103600 | 104700 |
| 7 | Devonport-Takapuna Local Board Area | 51600 | 52300 | 55700 | 58500 | 60500 | 62700 | 65700 | 68800 | 71800 | 74800 | 77600 |
| 8 | Henderson-Massey Local Board Area | 82200 | 90900 | 103700 | 113500 | 124600 | 134100 | 142500 | 150500 | 158200 | 165700 | 172900 |
| 9 | Waitākere Ranges Local Board Area | 42100 | 44200 | 47500 | 50700 | 54200 | 56900 | 58300 | 59500 | 60300 | 60900 | 61100 |
| 10 | Waitematā Local Board Area | 44900 | 52500 | 66600 | 81300 | 88500 | 94700 | 102800 | 111900 | 120700 | 129200 | 137400 |
| 11 | Whau Local Board Area | 59700 | 66200 | 73100 | 76700 | 84100 | 89500 | 95000 | 100500 | 105900 | 111000 | 115900 |
| 12 | Albert-Eden Local Board Area | 86100 | 90000 | 96200 | 100000 | 103700 | 106900 | 111200 | 116500 | 121600 | 126300 | 130600 |
| 13 | Puketāpapa Local Board Area | 46100 | 49800 | 53900 | 56300 | 60900 | 64600 | 68800 | 73000 | 77100 | 81000 | 84600 |
| 14 | Ōrākei Local Board Area | 70900 | 73000 | 78400 | 83700 | 87700 | 91300 | 96200 | 101100 | 105800 | 110300 | 114500 |
| 15 | Maungakiekie-Tāmaki Local Board Area | 61700 | 66400 | 70400 | 73700 | 80500 | 86600 | 93800 | 101300 | 108800 | 116100 | 123000 |
| 16 | Howick Local Board Area | 80900 | 97300 | 119100 | 135000 | 149400 | 160500 | 167800 | 174400 | 180400 | 185900 | 190900 |
| 17 | Māngere-Ōtāhuhu Local Board Area | 61400 | 64600 | 72500 | 75300 | 82700 | 88700 | 92500 | 96300 | 99800 | 103100 | 105900 |
| 18 | Ōtara-Papatoetoe Local Board Area | 67800 | 70400 | 76600 | 80300 | 90500 | 97700 | 101500 | 104000 | 106000 | 107500 | 108600 |
| 19 | Manurewa Local Board Area | 60100 | 69100 | 81600 | 87000 | 100900 | 109600 | 114000 | 118000 | 121700 | 124900 | 127600 |
| 20 | Papakura Local Board Area | 38000 | 38900 | 43100 | 48200 | 61100 | 70000 | 77000 | 81700 | 86100 | 90500 | 95000 |
| 21 | Franklin Local Board Area | 48400 | 52900 | 60600 | 68300 | 77700 | 85700 | 96900 | 109400 | 122000 | 134500 | 146900 |
d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()
(17800, 84495.2479338843, 190900)
%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
ims = []
for col in cols:
ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
ax.set_title(col)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("pop.gif", ims, fps=1)
with open('pop.gif','rb') as file:
display(Image(file.read()))
%%time
popcube = make_geocube(
vector_data=pop,
measurements=cols,
like=geocube, # Ensures dimensions match
fill=0 # NaNs, like offshore areas, will be 0
)
popcube
CPU times: user 1.98 s, sys: 65 µs, total: 1.98 s Wall time: 1.98 s
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
1996 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2001 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2006 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2013 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2018 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2023 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2028 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2033 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2038 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2043 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2048 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])for col in cols:
outfile = f"output/pop_{col}.tif"
if not os.path.isfile(outfile):
# byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
popcube[col].rio.to_raster(outfile, dtype=np.uint32)
%%time
roads = gpd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")
CPU times: user 8.58 s, sys: 9.64 ms, total: 8.59 s Wall time: 8.61 s
%%time
akl_roads = gpd.clip(roads, AKL)
CPU times: user 25.4 s, sys: 0 ns, total: 25.4 s Wall time: 25.4 s
# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway
| t50_fid | name_ascii | macronated | name | hway_num | rna_sufi | lane_count | way_count | status | surface | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 512 | 100120610 | KAIPARA COAST HIGHWAY | N | KAIPARA COAST HIGHWAY | 16 | 3007739 | 2 | None | None | sealed | LINESTRING (1732000.000 5944172.070, 1732048.5... |
| 2933 | 3198057 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748581.508 5968975.145, 1748558.4... |
| 2934 | 3198059 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748171.047 5971284.152, 1748129.9... |
| 3320 | 3200754 | PAERATA ROAD | N | PAERATA ROAD | 22 | 3000260 | 2 | None | None | sealed | LINESTRING (1767236.112 5888088.508, 1767244.3... |
| 3324 | 3200792 | UPPER HARBOUR MOTORWAY | N | UPPER HARBOUR MOTORWAY | 18 | 3047073 | 4 | None | None | sealed | LINESTRING (1747954.314 5927269.837, 1747970.0... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 138240 | 100048291 | AUCKLAND-WAIWERA MOTORWAY | N | AUCKLAND-WAIWERA MOTORWAY | 1 | 3067966 | 7 | None | None | sealed | LINESTRING (1755881.018 5922863.734, 1755886.4... |
| 138301 | 100048432 | AUCKLAND-HAMILTON MOTORWAY | N | AUCKLAND-HAMILTON MOTORWAY | 1 | 3017109 | 1 | None | None | sealed | LINESTRING (1765115.647 5909916.697, 1765092.7... |
| 138337 | 100048532 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 4 | None | None | sealed | LINESTRING (1748892.089 5949596.727, 1748892.0... |
| 138369 | 100048589 | PORT ALBERT ROAD | N | PORT ALBERT ROAD | 16 | 3013274 | 2 | None | None | sealed | LINESTRING (1734173.019 5980575.187, 1734175.6... |
| 138680 | 100118365 | SOUTH-WESTERN MOTORWAY | N | SOUTH-WESTERN MOTORWAY | 20 | 3018532 | 4 | None | None | sealed | LINESTRING (1760066.252 5908184.133, 1760043.9... |
426 rows × 11 columns
mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()
mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()
ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)
%%time
mway_cube = make_geocube(
vector_data=mway,
measurements=["lane_count"],
like=geocube, # Ensures dimensions match
fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube
CPU times: user 256 ms, sys: 157 µs, total: 256 ms Wall time: 254 ms
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
lane_count (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
(1320, 1001)
Text(0.5, 1.0, 'Distance (m)')
%%time
skytower = [-36.84838748948485, 174.7621736911587]
skytower = gpd.points_from_xy(x=[skytower[1]], y=[skytower[0]])
cbd = gpd.GeoDataFrame([{"name": "Skytower", "value": 1}], geometry=skytower, crs="EPSG:4326").to_crs(epsg=2193)
display(cbd)
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
vector_data=cbd,
like=geocube, # Ensures dimensions match
fill=0
)
display(cbd_cube)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
cbd_cube.value.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
| name | value | geometry | |
|---|---|---|---|
| 0 | Skytower | 1 | POINT (1757109.057 5920497.145) |
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
value (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])(1320, 1001) CPU times: user 535 ms, sys: 408 µs, total: 536 ms Wall time: 535 ms
Text(0.5, 1.0, 'Distance (m)')
bounds = AKL.total_bounds.tolist()
bounds
[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]
zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles
['EJ.tif', 'DM.tif', 'EL.tif', 'DL.tif', 'DJ.tif', 'FK.tif', 'DK.tif', 'EK.tif', 'FL.tif']
tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)
print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')
-32767 -18316.517014943143 697 (1, 1320, 1018)
<AxesSubplot:>
cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape
(1, 1320, 1001)
NODATA = -32767
warped_DEM = np.full_like(cbd_dist.read(), NODATA, dtype=np.int16)
rasterio.warp.reproject(DEM,
warped_DEM,
src_transform=transform,
src_nodata=NODATA,
dst_nodata=NODATA,
src_crs=tile_datasets[0].crs,
dst_crs=cbd_dist.crs,
dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)
(1, 1320, 1018) (1, 1320, 1001)
meta = tile_datasets[0].meta
print(meta)
meta.update({
"dtype": "int16",
"height": cbd_dist.height,
"width": cbd_dist.width,
"transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
with rasterio.open(outfile, "w", **meta) as dest:
dest.write(warped_DEM)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
0.0, -100.0, 6002400.0)}
gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(np.where(slope>=0, slope, np.nan), cmap='terrain')
-9999.0 -5778.3354 52.492313 (1, 1320, 1001)
<AxesSubplot:>
rasters = [rasterio.open(f"output/{f}") for f in os.listdir("output") if f.endswith(".tif")]
params = set([(r.shape, r.res, r.crs, r.count, r.bounds) for r in rasters])
print(params)
# Assert all rasters have the same shape, pixel size, CRS, number of bands, and bounds
assert len(params) == 1
{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}
files = []
for file in sorted(os.scandir("output"), key=lambda file: file.name):
files.append({
"name": file.name,
"filesize (MB)": round(file.stat().st_size / 1024 / 1024, 2),
"last modified": datetime.datetime.fromtimestamp(file.stat().st_mtime)
})
output = pd.DataFrame(files)
display(output)
print(f"Total: {output['filesize (MB)'].sum().round()}MB")
| name | filesize (MB) | last modified | |
|---|---|---|---|
| 0 | cbd.tif | 1.26 | 2021-05-04 14:19:19.670 |
| 1 | cbd_dist.tif | 2.52 | 2021-05-04 14:19:19.720 |
| 2 | dem.tif | 2.52 | 2021-05-04 14:19:22.890 |
| 3 | land_use_1996.tif | 1.26 | 2021-05-04 14:45:25.500 |
| 4 | land_use_2001.tif | 1.26 | 2021-05-04 14:45:25.500 |
| 5 | land_use_2008.tif | 1.26 | 2021-05-04 14:45:25.510 |
| 6 | land_use_2012.tif | 1.26 | 2021-05-04 14:45:25.510 |
| 7 | land_use_2018.tif | 1.26 | 2021-05-04 14:45:25.520 |
| 8 | mway.tif | 1.26 | 2021-05-04 14:19:18.360 |
| 9 | mway_dist.tif | 2.52 | 2021-05-04 14:19:18.900 |
| 10 | pop_1996.tif | 5.04 | 2021-05-04 14:18:13.670 |
| 11 | pop_2001.tif | 5.04 | 2021-05-04 14:18:13.680 |
| 12 | pop_2006.tif | 5.04 | 2021-05-04 14:18:13.690 |
| 13 | pop_2013.tif | 5.04 | 2021-05-04 14:18:13.700 |
| 14 | pop_2018.tif | 5.04 | 2021-05-04 14:18:13.710 |
| 15 | pop_2023.tif | 5.04 | 2021-05-04 14:18:13.720 |
| 16 | pop_2028.tif | 5.04 | 2021-05-04 14:18:13.730 |
| 17 | pop_2033.tif | 5.04 | 2021-05-04 14:18:13.740 |
| 18 | pop_2038.tif | 5.04 | 2021-05-04 14:18:13.740 |
| 19 | pop_2043.tif | 5.04 | 2021-05-04 14:18:13.760 |
| 20 | pop_2048.tif | 5.04 | 2021-05-04 14:18:13.770 |
| 21 | slope.tif | 5.04 | 2021-05-04 14:19:22.930 |
Total: 77.0MB